Detection of synchronization from univariate data using wavelet transform 
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A method is proposed for detecting from univariate data the presence of synchronization of a 
self-sustained oscillator by external driving with varying frequency. The method is based on the 
analysis of difference between the oscillator instantaneous phases calculated using continuous wavelet 
transform at time moments shifted by a certain constant value relative to each other. We apply our 
method to a driven asymmetric van der Pol oscillator, experimental data from a driven electronic 
oscillator with delayed feedback and human heartbeat time series. In the latest case, the analysis 
of the heart rate variability data reveals synchronous regimes between the respiration and slow 
oscillations in blood pressure. 
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I. INTRODUCTION 

Detecting regimes of synchronization between self- 
sustained oscillators is a typical problem in studying their 
interaction. Two types of interaction are generally rec- 
ognized P, H, S 01- The first one is a unidirectional 
coupling of oscillators. It can result in synchronization 
of a self-sustained oscillator by an external force. In this 
case the dynamics of the oscillator generating the driv- 
ing signal does not depend on the driven system behav- 
ior. The second type is a mutual coupling of oscillators. 
In this case the interaction can be more effective in one 
of the directions, approaching in the limit to the first 
type, or can be equally effective in both directions. In 
the event of mutual coupling, synchronization is the re- 
sult of the adjustment of rhythms of interacting systems. 
To detect synchronization one can analyze the ratio of 
instantaneous frequencies of interacting oscillators and 
the dynamics of the generalized phase difference [3| . As 
a quantitative characteristic of synchronization one can 
use the phase synchronization index d 0] or the measure 
of synchronization 0, [H . 

Synchronization of interacting systems including the 
chaotic ones has been intensively studied in recent years. 
The main ideas in this area have been introduced using 
standard models 0, 0, i i, 0, S H M M El El E3 ) • 
At present, more attention is focused on application of 
the developed techniques to living systems. In particular, 
much consideration is being given to investigation of syn- 



chronization between different brain areas [a. Il5l. Il6l. Il7| 
and to studying synchronization in the human cardiores- 
piratory system [18|, [11 [H [22j ■ Investigating such 



systems one usually deals with the analysis of short time 
series heavily corrupted by noise. In the presence of noise 
it is often difficult to detect the transitions between syn- 
chronous and nonsynchronous regimes. Besides, even in 
the region of synchronization a 27r-phase jumps in the 
temporal behavior of the generalized phase difference can 
take place. Moreover, the interacting systems can have a 
set of natural rhythms. That is why it is desirable to an- 
alyze synchronization and phase locking at different time 

scales m M M M El. 
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A striking example of interaction between various 
rhythms is the operation of the human cardiovascular 
system (CVS). The main rhythmic processes governing 
the cardiovascular dynamics are the main heart rhythm, 
respiration, and the process of slow regulation of blood 
pressure and heart rate havin g in humans the fundamen- 
tal frequency close to 0.1 Hz [26]. Owing to interaction, 
these rhythms appear in various signals: electrocardio- 
gram (ECG), blood pressure, blood flow, and heart rate 
variability (HRV) [27j. Recently, it has been found that 
the main rhythmic processes operating within the CVS 
can be synchronized [H El HI HH • It has been shown 
that the systems generating the main heart rhythm and 
the rhythm associated with slow oscillations in blood 
pressure can be regarded as self-sustained oscillators, and 
that the respiration can be regarded as an external forc- 
ing of these systems [H [H| ■ 

Recently, we have proposed a method for detecting the 
presence of synchronization of a self-sustained oscillator 
by external driving with linearly varying frequency [22| . 
This method was based on a continuous wavelet trans- 
form of both the signals of the self-sustained oscillator 
and external force. However, in many applications the 
diagnostics of synchronization from the analysis of uni- 
variate data is a more attractive problem than the de- 
tection of synchronization from multivariate data. For 
instance, the record of only a univariate signal may be 
available for the analysis or simultaneous registration of 
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different variables may be rather difficult. In this paper 
we propose a method for detection of synchronization 
from univariate data. However, a necessary condition 
for application of our method is the presence of a driv- 
ing signal with varying frequency. For the mentioned 
above cardiovascular system our method gives a possibil- 
ity to detect synchronization between its main rhythmic 
processes from the analysis of the single heartbeat time 
series recorded under paced respiration. 

The paper is organized as follows. In Sec. II we de- 
scribe the method for detecting synchronization from uni- 
variate data. In Sec. Ill the method is tested by applying 
it to numerical data produced by a driven asymmetric 
van der Pol oscillator. In Sec. IV the method is used 
for detecting synchronization from experimental time se- 
ries gained from a driven electronic oscillator with de- 
layed feedback. Section V presents the results of the 
method application to studying synchronization between 
the rhythms of the cardiovascular system from the anal- 
ysis of the human heart rate variability data. In Sec. VI 
we summarize our results. 



II. METHOD DESCRIPTION 

Let us consider a self-sustained oscillator driven by ex- 
ternal force T with varying frequency 



x = H(x) + eF($(t)), 



(1) 



where H is the operator of evolution, e is the driving 
amplitude, and <f>(i) is the phase of the external force 
defining the law of the driving frequency u>d(t) variation: 



dt 



(2) 



In the simplest case the external force is described by a 
harmonic function T(Q(t)) = sin$(i). 

Assume that we have at the disposal a univariate time 
series x{t) characterizing the response of the oscillator 
to the driving force T . Let us define from this time series 
the phase ipo (t) of oscillations at the system (p} basic fre- 
quency /o- The main idea of our approach for detecting 
synchronization from univariate data is to consider the 
temporal behavior of the difference between the oscillator 
instantaneous phases at the time moments t and t + r. 
We calculate the phase difference 



A(p (t) = <p a (t + t) - (fo{t), 



(3) 



where r is the time shift that can be varied in a wide 
range. Note, that <po(t) and (po(t + r) are the phases 
of the driven self-sustained oscillator corresponding to 
oscillations at the first harmonic of the oscillator basic 
frequency / . 

The variation of driving frequency is crucial for the 
proposed method. Varying in time, the frequency of 
the external force sequentially passes through the regions 
of synchronization of different orders 1:1,2:1,..., 



n : 1, n : m, . ..(n,m = 1,2,3,...). Within the 
time intervals corresponding to asynchronous dynamics 
the external signal practically has no influence on the 
dynamics of the basic frequency fo in the oscillator (Q]) 
spectrum. Thus, the phase of oscillator varies linearly 
outside the regions of synchronization, (fo(t) — 2irfot + (p, 
where (p is the initial phase. Then, from Eq. ([3]) it follows 



(4) 



i.e., the phase difference Aipo(t) is constant within the 
regions of asynchronous dynamics. 

Another situation is observed in the vicinity of the 
time moments t ns where the driving frequency uj^it) « 
(2irn/m)fo and n : m synchronization takes place. For 
simplicity let us consider the case of 1 : 1 synchronization. 
In the synchronization (Arnold) tongue the frequency of 
the system ([T]) nonautonomous oscillations is equal to the 
frequency ([2]) of the external force and the phase differ- 
ence between the phase of the driven oscillator (po (t) and 
the phase &(t) of the external force, A<j)(t) = ip (t) — $(t), 
is governed in a first approximation by the Adler equa- 
tion [28|. It follows from the Adler equation that in the 
region of 1 : 1 synchronization the phase difference A(f>(t) 
varies by n. 

Representing the driven oscillator phase as <po(t) = 
A<j)(t) + $(t), we obtain from Eq. (J3J): 



A<p (t) =$(< + t) -$(*)+ 7, 



(5) 



where 7 = A(fr(t + r) — A(f>(t) s» const is the correction of 
the phase difference that appears due to synchronization 
of the system by external force. Expanding the phase 
$(f + t) in a Taylor series we obtain 



Aifi (t) = 7 + —^t + -^A±t 2 



dt ' ' 2 dt 2 
Taking into account Eq. ([2]) we can rewrite Eq. ([6]) as 



(0) 



Aif Q {i) = 7 + LO d {t)T + — 



(J) 



Thus, the behavior of the phase difference ([3]) is defined 
by the law of the driving frequency u>d(t) variation. 

For the linear variation of the driving frequency, 
oj d (t) =a + fit, from Eq. © it follows 



A(f (t) = 7 + ar + (3t 2 /2 + rf3t. 



(8) 



Consequently, in the region of synchronization the phase 
difference varies linearly in time, Atpo(t) ~ t. In the 
case of the nonlinear variation of ujd(t), the dynamics 
of Aipoit) is more complicated. However, if u>d(t) varies 
in a monotone way and the time of its passing through 
the synchronization tongue is small, one can neglect the 
high-order terms of the expansion and consider the law 
of Aipo(t) variation as the linear one. We will show below 
that this assumption holds true for many applications. 
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The absolute value of the change in the phase differ- 
ence Aipo(t) within the synchronization region can be 
estimated using Eq. J7|: 

Aip s = A(p (t 2 ) - Atp (ti) = (ujdfa) - LOd(ti))r+ 



duj d (t) 



dt 



du> d (t) 



dt 




(9) 



where t\ and ti are the time moments when the fre- 
quency of the external force passes through, respectively, 
the low-frequency and high-frequency boundaries of the 
synchronization tongue. Assuming that the rate of u) d {t) 
variation is slow, we can neglect the terms containing the 
derivatives of u>d(t) and obtain 



Aip s s» Awt, 



(10) 



where Auj = ^(^2) ~^>d(ti) is the bandwidth of synchro- 
nization. 

The obtained estimation corresponds to the case of 
1 : 1 synchronization, characterized by equal values of 
the driving frequency fd and the oscillator frequency /o, 
fd/fo — 1< However, the considered approach can be 
easily extended to a more complicated case of n : m syn- 
chronization. In this case the change in Aipa(t) within 
the region of synchronization takes the value 



Tfi 

Atp s = — Aujt. 
n 



(11) 



Hence, the analysis of the phase difference © behav- 
ior allows one to distinguish between the regimes of syn- 
chronous and asynchronous dynamics of driven oscillator. 
The phase difference A<po(t) is constant for the regions 
of asynchronous dynamics and demonstrates monotone 
(often almost linear) variation by the value Aip s defined 
by Eq. (|lip within the regions of synchronization. 

To define the phase (po(t) of oscillations at the basic 
frequency we use the approach based on the continuous 
wavelet transform 0, [a, [241 |29| . It is significant, that 
the wavelet transform [30, [31| is the powerful tool for the 
analysis of nonlinear dynamical system behavior. The 
continuous wavelet analysis has been applied in the stud- 
ies of phase synchronization of chaotic neural oscillations 
in the brain [32], [H, [34], [H, Hi| , electroencephalogram sig- 
nals [37j , R-R intervals and arterial blood pressure oscil- 
lations in brain injury (38j . chaotic laser array (39| . It has 
also been used to detect the main frequency of the oscil- 
lations in nephron autoregulation (40| and coherence be- 
tween blood flow and skin temperature oscillations [4lj . 
In these recent studies a continuous wavelet transform 
with various mother wavelet functions has been used for 
introducing the instantaneo us p hases of analyzed signals. 
In particular, in Refs. [34], l37f ] a comparison of Hilbcrt 
transform and wavelet method with the mother Morlet 
wavelet has been carried out and good conformity be- 
tween these two methods has been shown for the analysis 
of neuronal activity. It is important to note, that in all 



the above mentioned studies the wavelet transform has 
been used for the analysis of synchronization from bivari- 
ate data, when the generalized phase difference Aip(t) of 
both analyzed rhythms was investigated. The proposed 
method allows one to detect synchronization from the 
analysis of only the one signal of the oscillator response to 
the external force with monotonically varying frequency. 
Taking into account the high efficiency of the analysis of 
synchronization with the help of the continuous wavelet 
transform using bivariate data, we will use the continu- 
ous wavelet transform for determining the instantaneous 
phase of the analyzed univariate signal. 

The continuous wavelet transform [3(| HH of the signal 
x{t) is defined as 



w( S ,t )= / X (t)r s , t0 (t)dt, 



(12) 



where ip a ,to W i s t ne wavelet function related to the 
mother wavelet ipo(t) as ip s ,t (t) — (lAA) V'o ((t — t )/s). 
The time scale s corresponds to the width of the wavelet 
function, to is the shift of the wavelet along the time axis, 
and the asterisk denotes complex conjugation. It should 
be noted that the wavelet analysis operates usually with 
the time scale s instead of the frequency /, or the cor- 
responding period T = 1//, traditional for the Fourier 
transform. 

The wavelet spectrum 



W(s,t ) = \W(s,t )\exp\j<p s {t )] 



(13) 



describes the system dynamics for every time scale s at 
any time moment t$. The value of |W(s,io)| determines 
the presence and intensity of the time scale s at the 
time moment to- We use the complex Morlet wavelet 
[42l ] V'o(^) = exp[jar]] exp [— rj 2 /2] as the mother 

wavelet function. The choice of the wavelet parameter 
a = 2n provides the simple relation / sa 1/s between the 
frequency / of the Fourier transform and the time scale 

s M- 



III. METHOD APPLICATION TO DETECTING 
SYNCHRONIZATION IN A DRIVEN 
ASYMMETRIC VAN DER POL OSCILLATOR 

A. Model 

Let us consider the asymmetric van der Pol oscillator 
under external force with linearly increasing frequency: 



(l — fix 



x z )x + 



n 2 x = £sin$(i), (14) 



where /1 is the parameter characterizing the system asym- 
metry, = 0.247T is the natural frequency, and e and <&(i) 
are, respectively, the amplitude and phase of the exter- 
nal force. The phase = 2tt [(a + (3t/T)} t defines the 
linear dependence of the driving frequency Lo d (t) on time: 



Ud(t) 



d${t) 
dt 



2ir [a + 2f3t/T] 



(15) 



4 



where a — 0.03, /3 = 0.17, and T — 1800 is the maximal 
time of computation. This system has been considered 
in Ref. [13] as a model for studying synchronization be- 
tween the respiration, which can be regarded as an ex- 
ternal force, and the process of slow regulation of blood 
pressure and heart rate, which can be treated as a self- 
sustained oscillator. In the present paper we use this 
model system for testing our new method of detecting 
synchronization from univariate data. The chosen values 
of the model parameters provide close correspondence of 
frequencies and the ways of the driving frequency varia- 
tion in the simulation and experimental study described 
in Sec. V. The parameter /i is chosen to be equal to unity 
throughout this paper. In this case the phase portrait of 
oscillations is asymmetric and the power spectrum con- 
tains both odd and even harmonics of the basic frequency 
fo = 0.0973, as well as the power spectrum of the low- 
frequency fluctuations of blood pressure and heart rate 
[22I ] . Recall that the classical van der Pol oscillator with 
fi = has a symmetric phase portrait and its power spec- 
trum exhibits only odd harmonics of fo- We calculate the 
time series of nonautonomous asymmetric van der Pol os- 
cillator (|14p at e — 0.2 using a fourth-order Runge-Kutta 
method with the integration step At = 0.01. 



B. Results 

Fig. Q] shows the amplitude spectrum |VF(s,io)| 01 the 
wavelet transform for the signal of driven oscillator (fl~4|) . 
The Morlet wavelet is used as the mother wavelet func- 
tion throughout the paper. The wavelet parameter is 
chosen to be a = 2ir, unless otherwise specified. The 
time scale sq corresponding to the first harmonic of the 
oscillator basic frequency fo is indicated in Fig. [1] by the 
dot-and-dash line. The dashed line indicates the time 
scale Si corresponding to the linearly increasing driving 
frequency u>d(t). The analysis of the wavelet power spec- 
trum reveals the classical picture of oscillator frequency 
locking by the external driving. As the result of this 
locking, the breaks appear close to the time moments t s 
and t2s denoted by arrows, when the driving frequency 
is close to the oscillator basic frequency (uJd{t s ) « 2nfo) 
or to its second harmonic (tudfas) ~ 47r/o), respectively. 
These breaks represent the entrainment of oscillator fre- 
quency and its harmonic by external driving. If the de- 
tuning S = (ud — 27r/o) is great enough, the frequency of 
oscillations returns to the oscillator basic frequency. 

The dynamics of the phase differences Aipo(t) deter- 
mined by Eq. (J3J) is presented in Fig. [5^ for different 
positive r values. One can see in the figure the regions 
where A(po(t) is almost constant. These are the regions 
of asynchronous dynamics, when the driving frequency is 
far from the oscillator basic frequency and its harmon- 
ics. The regions of monotone increase of A(fo(t) are also 
well-pronounced in Fig. [2^. These are the regions of syn- 
chronization observed in the vicinity of the time moments 
t ns , when ajd(t ns ) ~ lim fo- 
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FIG. 1: (Color online) Shaded plot of the wavelet power spec- 
trum IVF^to)! for the signal generated by oscillator (|14|l. 
Time is shown on the abscissa and time scale is shown on the 
ordinate. The color intensity is proportional to the absolute 
value of the wavelet transform coefficients. The values of the 
coefficients are indicated by the scale from the right side of 
the figure. 
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FIG. 2: (Color online) Phase differences Atpo(t) Q calculated 
at the time scale so corresponding to the basic frequency fo = 
0.0973 of the driven asymmetric van der Pol oscillator (|14|l for 
different r > (a) and r < (b). 



The proposed method offers several advantages over 
the method in Ref. [22| based on the analysis of the phase 
difference between the signals of oscillator and the exter- 
nal force. First, the regions of Atpo(t) monotone variation 
corresponding to synchronous regimes are easily distin- 
guished from the regions of constant Aipo(t) value cor- 
responding to asynchronous dynamics. Second, the new 
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method is considerably more sensitive than the previous 
one because the phase difference is examined at the time 
scales having high amplitude in the wavelet spectrum. In 
particular, the region of 3 : 1 synchronization in the vicin- 
ity of the time moment t^ s denoted by arrow is clearly 
identified in Fig. O Third, the proposed method is sub- 
stantially simpler than the method of the phase difference 
calculation along the scale varying in time [13] ■ 

It follows from Eq. ([7]) that in the region of synchro- 
nization the change of the phase difference Atp (t) in- 
creases with r increasing. As the result, the presence 
of interval of Aipo(t) monotone variation becomes more 
pronounced, Fig. [2^i. This feature helps to detect the ex- 
istence of synchronization especially in the case of high- 
order synchronization and noise presence. However, the 
accuracy of determining the boundaries of the region of 
synchronization decreases as r increases. 

It should be noted that for negative r values the mono- 
tone reduction of the phase difference is observed in the 
region of synchronization, Fig. [3b. As it can be seen 
from Fig. [2b, the increase of r by absolute value leads to 
increase of A^o (t) variation in the region of synchroniza- 
tion as well as in the case of positive r. 



C. Influence of noise and inaccuracy of the basic 
time scale definition 

Experimental data, especially those obtained from liv- 
ing systems, are always corrupted by noise. Besides, in 
many cases it is not possible to define accurately the basic 
frequency of the system under investigation. For exam- 
ple, interaction between the human cardiovascular and 
respiratory systems and nonstationarity hampers accu- 
rate estimation of natural frequencies for cardiovascular 
rhythms. Therefore, the actual problem is to test the 
method efficiency for detecting synchronization in the 
presence of additive noise and inaccuracy of the basic 
frequencies estimation. 

To analyze the influence of noise on the diagnostics of 
synchronization we consider the signal 



x n (t) = x(t) + DC(t), 



(16) 



where x(t) is the signal of the asymmetric van der Pol 
oscillator ([H)) , £(£) is the additive noise with zero mean 
and uniform distribution in the interval [—0.5, 0.5], and 
D is the intensity of noise. To simulate the noisy signal 
£(f) we use the random- number generator described in 
Ref. [H. 

Typical time series x n it) generated by Eq. ([To]) for dif- 
ferent intensities of noise are presented in Fig. [3ji for the 
region of 1 : 1 synchronization. In spite of the signifi- 
cant distortion of the signal by noise its wavelet power 
spectrum, Fig. [3b, still allows to reveal the main features 
of the system dynamics. In particular, the dynamics of 
the time scale so and the effect of frequency entrainment 
in the region of 1 : 1 synchronization indicated by arrow 
are recognized in Fig. [3b- Hence, the use of the wavelet 
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FIG. 3: (Color online) (a) Parts of the time series of the 
signal (|16[) for different intensities D of additive noise, (b) 
Wavelet power spectrum |W(s, to) | of the signal x n (t) at the 
noise intensity D = 10. The dot-and-dash line indicates the 
time scale so corresponding to the oscillator basic frequency 
fo- (c, d) Phase differences Aipo{t) for different intensities D 
of noise at r = 10 (c) and r = 100 (d). The inset in (c) is the 
enlarged fragment of the region of 1 : 1 synchronization. 



transform for determining the phases of the signal and its 
harmonics allows one to detect the regimes of synchro- 
nization from noisy time series. 

The phase differences Aip (t) calculated using Eq. ([3J 
with t = 10 are shown on Fig. [3J: for different inten- 
sities D of additive noise. The dependence Aipo(t) be- 
comes more jagged as D increases. However, for D < 10 
we can identify the regions where the phase difference 
demonstrates near-monotone variation. In the average 
this variation is about the same as in the case of noise 
absence (see the inset in Fig. [3b). Fig. [3bl shows Aip {t) 
for r = 100. In this case it is possible to detect the pres- 
ence of synchronization for significantly higher levels of 
noise than in the case of small r. The reason is that the 
value of Aip s (jlip increases in the region of synchroniza- 
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Atpo [radj 




FIG. 4: (Color online) Phase differences Atpo(t) calculated 
at the time scales Si = so + As for r = 100 and D — 10. 
The curve numbers correspond to the following time scaled: 
(1) si = 7.28 < s , (2) si = 8.28 < s , (3) si = s = 10.28, 
(4) Sl = 12.28 > s , (5) si = 15.28 > s . 



tion as the time shift r increases, whereas the amplitude 
of Aipo(t) fluctuations caused by noise does not depend 
on r. For very large intensities of noise (D = 50 in Fig. [3]) 
the synchronous behavior is not so clearly pronounced as 
at smaller D values, but it should be mentioned that in 
this case the power of noise exceeds the power of the 
oscillator signal in several times. 

Let us consider the method efficiency in the case when 
the scale s of observation differs from the time scale so 
associated with the oscillator basic frequency Jq. Fig. Q] 
illustrates the behavior of the phase difference Aipo(t) 
calculated for the time series of Eq. flU]) at the time 
scales si = s + As, where As is the detuning of the 
scale with respect to the basic scale so ~ l//o = 10.28. 
It can be seen from the figure that for |As| < 2.0 the 
phase dynamics is qualitatively similar to the case of ac- 
curate adjustment of the scale s to the basic scale So- 
At greater As values the phase difference demonstrates 
significant fluctuations impeding to detect the epochs of 
A<po(t) monotone variation. Thus, to detect synchroniza- 
tion using the proposed method one needs to know only 
approximately the basic time scale sq. 



IV. INVESTIGATION OF SYNCHRONIZATION 
IN A DRIVEN ELECTRONIC OSCILLATOR 
WITH DELAYED FEEDBACK 

A. Experiment description 

We apply the method to experimental data gained 
from a driven electronic oscillator with delayed feedback. 
A block diagram of the experimental setup is shown in 



Fig. [5l The oscillator represents the ring system com- 
posed of nonlinear, delay, and inertial elements. The 
role of nonlinear clement is played by an amplifier with 
the quadratic transfer function. This nonlinear device is 
constructed using bipolar transistors. The delay line is 
constructed using digital elements. The inertial proper- 
ties of the oscillator are defined by a low-frequency first- 
order i?C-filter. The analogue and digital elements of the 
scheme are connected with the help of analog-to-digital 
(ADC) and digital-to-analog converters (DAC). To gen- 
erate the driving signal we use the sine- wave generator 2 
whose frequency is modulated through the wobble input 
by the signal of the sawtooth pulse generator 1. The 
driving signal is applied to the oscillator using the sum- 
mator E. The considered oscillator is governed by the 
first-order time-delay differential equation 

RCU(t) = -U{t) + F(U(t - d)) + U sin(2TT f ext (t)t), 

where U(t) and U(t — d) are the delay line input and 
output voltages, respectively, d is the delay time, R and C 
are the resistance and capacitance of the filter elements, 
F is the transfer function of the nonlinear device, Uq 
is the amplitude of the driving signal, and f ex t is the 
driving frequency. We record the signal U(t) using an 
analog-to-digital converter with the sampling frequency 
/ = 15 kHz at d = 1.5 ms and RC — 0.46 ms under the 
following variation of the driving frequency 

f ex t(t)=vl0 u ^ /2 , (18) 

where v = 220 Hz and the control voltage U w (t) varies 
linearly from V to 16 V within 800 ms providing f ext 
variation from 220 Hz to 1000 Hz. Under the chosen pa- 
rameters the considered oscillator demonstrates periodic 
oscillations with the period T = 3.7ms. Four exper- 
iments were carried out at different amplitudes of the 
external driving equal to 0.5 V, IV, 1.5 V, and 2 V. The 
amplitude of driven oscillation was about 3 V. 

B. Results 

The experimental time series of the electronic oscillator 
with delayed feedback driven by the external force with 
varying frequency (fl"5)) are depicted in Fig. \E\ for two val- 
ues of the driving amplitude. The results of investigation 
of the oscillator synchronization by the external driving 
are presented in Fig. [7] The phase differences Atp (t) 
defined by Eq. ([3]) are calculated under different driv- 
ing amplitudes Uq for the time shift r = —0.66 ms. One 
can clearly identify in the figure the regions of Atpo(t) 
monotone variation corresponding to the closeness of the 
driving frequency to the oscillator basic frequency and its 
harmonics. These regions of synchronous dynamics are 
indicated by arrows. 

It is well seen from Fig.[7]that the interval of monotone 
variation of Aipo(t) increases with increasing amplitude 
of the driving force. This fact agrees well with the known 
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FIG. 5: Block diagram of the electronic oscillator with delayed feedback driven by the signal with varying frequency. 
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FIG. 6: (Color online) Time series of electronic oscillator 
with delayed feedback under external driving with varying 
frequency (|18p and the driving amplitude Uo = 0.5 V (a) and 
C/ = 2V(b). 
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FIG. 7: (Color online) Phase differences A(po(t) © calcu- 
lated at the time scale so corresponding to the basic frequency 
fo = 270 Hz of the driven electronic oscillator with delayed 
feedback. The curve numbers correspond to different ampli- 
tudes Uo of the external force: (1) U = 0.5 V, (2) U = 1 V, 
(3) Uo = 1.5 V, (4) Uo = 2V. 



V. SYNCHRONIZATION OF SLOW 
OSCILLATIONS IN BLOOD PRESSURE BY 
RESPIRATION FROM THE DATA OF HEART 
RATE VARIABILITY 



effect of extension of the region of synchronization with 
increase in the amplitude of the external driving. Note, 
that in spite of the nonlinear variation of the driving fre- 
quency, at small driving amplitudes the phase difference 
A<po(t) varies almost linearly in time within the synchro- 
nization tongue as it was discussed in Sec. [ill For the 
large driving amplitude (Uo = 2 V) the synchronization 
tongue is wide enough and the phase difference behavior 
begins to depart from linearity. Nevertheless, the varia- 
tion of Atpo(t) remains the monotone one and allows us 
to detect the presence of synchronization and estimate 
the boundaries of the synchronization tongue. 



In this section we investigate synchronization between 
the respiration and rhythmic process of slow regulation 
of blood pressure and heart rate from the analysis of uni- 
variate data in the form of the heartbeat time series. This 
kind of synchronization has been experimentally studied 
in [2l|, [z| 0, EH • We studied eight healthy volunteers. 
The signal of ECG was recorded with the sampling fre- 
quency 250 Hz and 16-bit resolution. Note, that accord- 
ing to [40| the sampling frequency 250 Hz used in our 
experiments suffices to detect accurately the time mo- 
ment of R peak appearance. The experiments were car- 
ried out under paced respiration with the breathing fre- 
quency linearly increasing from 0.05 Hz to 0.3 Hz within 
30 min. We specially included the lower frequencies for 
paced respiration in order to illustrate the presence of the 
most pronounced regime of 1:1 synchronization between 
the respiration and slow oscillations in blood pressure. 
The rate of respiration was set by sound pulses. The de- 
tailed description of the experiment is given in Ref. [2l| . 
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Extracting from the ECG signal the sequence of R- 
R intervals, i.e., the series of the time intervals between 
the two successive R peaks, we obtain the information 
about the heart rate variability. The proposed method of 
detecting synchronization from uniform data was applied 
to the sequences of R-R intervals. 

A typical time series of R-R intervals for breathing at 
linearly increasing frequency is shown in Fig. [5^,. Since 
the sequence of R-R intervals is not equidistant, we ex- 
ploit the technique for applying the continuous wavelet 
transform to nonequidistant data. The wavelet spectra 
| W(s, *o) | for different parameters a of the Morlet wavelet 
are shown in Figs. [SJd and [5J; for the sequence of R-R 
intervals presented in Fig. [SK For greater a values the 
wavelet transform provides higher resolution of frequency 
[3ll | and better identification of the dynamics at the time 
scales corresponding to the basic frequency of oscillations 
and the varying respiratory frequency. In the case of 
a = 2tt the time scale s of the wavelet transform is very 
close to the period T of the Fourier transform and the 
values of s are given in seconds in Fig.[5jD. Generally, the 
time scale s is related to the frequency / of the Fourier 
transform by the following equation: 



R-R [s] 

0.9 



a + Vtr 2 + 2 

4^7 • 



(19) 



Because of this, the units on the ordinates are different 
in Figs. [5Jd and[8j:. The wavelet spectra in these figures 
demonstrate the high-amplitude component correspond- 
ing to the varying respiratory frequency manifesting itself 
in the HRV data. The self-sustained slow oscillations in 
blood pressure (Mayer wave) have in humans the basic 
frequency of about 0.1 Hz, or respectively, the basic pe- 
riod close to 10 s. The power of this rhythm in the HRV 
data is less than the power of respiratory oscillations. As 
the result, the time scale sq is weakly pronounced in the 
spectra. 

Fig. [9] presents the phase differences Aipo(t) calculated 
for R-R intervals of four subjects under respiration with 
linearly increasing frequency. All the curves in the figure 
exhibit the regions with almost linear in the average vari- 
ation of Atpo(t) indicating the presence of synchronous 
dynamics. In particular, the region of 1 : 1 synchroniza- 
tion is observed within the interval 200-600 s when the 
frequency of respiration is close to the basic frequency of 
the Mayer wave. This region is marked by arrow. In this 
region the frequency of blood pressure slow oscillations 
is locked by the increasing frequency of respiration and 
increases from 0.07 Hz to 0.14 Hz. Outside the interval of 
synchronization, t < 200 s and t > 600 s, the phase differ- 
ences demonstrate fluctuations caused by the high level 
of noise and nonstationarity of the experimental data. 
Some of these fluctuations take place around an average 
value as well as in the case of the driven van der Pol os- 
cillator affected by noise (see Fig. [3]). The frequency of 
blood pressure slow oscillations demonstrates small fluc- 
tuations around the mean value of about 0.1 Hz outside 
the interval of synchronization. 
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FIG. 8: (Color online) Typical sequence of R-R intervals for 
the case of breathing with linearly increasing frequency (a) 
and its wavelet power spectra at a — 2ir (b) and a = 16 (c). 
The dashed lines indicate the time scale so corresponding to 
the basic frequency /o = 0.1 Hz of slow oscillations in blood 
pressure. 



The phase differences in Fig. [5^, are plotted for differ- 
ent t. As the time shift r increases, so does the range of 
Aipo(t) monotone variation in the region of synchroniza- 
tion. This result agrees well with the results presented in 
Sec. III. Similar behavior of Aipo(t) is observed for each of 
the eight subjects studied. In Fig. E^b) phase differences 
AtfQ (t) computed for R-R intervals of another three sub- 
jects are presented. The phase differences demonstrate 
the wide regions of almost linear variation for all the sub- 
jects. Such behavior of the considered phase difference 
cannot be observed in the absence of synchronization, 
if only the modulation of blood pressure oscillations by 
respiration is present. These results allow us to confirm 
the conclusion that the slow oscillations in blood pressure 
can be synchronized by respiration. However, to come to 
this conclusion, the proposed method needs only univari- 
ate data in distinction to the methods [2l|, [22j based on 
the analysis of bivariate data. Note, that paper [2l[ con- 
tains the more detailed investigation of synchronization 
between the respiration and slow oscillations in blood 
pressure than the present one. Recent reports (see, for 
examples, [13, IH, [5(3] ) focused on examining the relation- 
ship between respiration and heart rate have shown that 
there is nonlinear coupling between respiration and heart 
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FIG. 9: (Color online) Phase differences Ayo(i) calculated 
at the time scale so corresponding to the basic frequency 
fo = 0.1 Hz of the Mayer wave, (a) Phase differences com- 
puted at different time shifts r for R-R intervals of one of 
the subjects. The curve numbers correspond to different time 
shifts: (1) r = 30 s, (2) r = 50 s, (3) r = 100 s. (b) Phase 
differences computed for R-R intervals of the other three sub- 
jects. 

rate. In particular, such coupling is well studied for the 
respiratory modulation of heart rate [i^, [5(| known as 
respiratory sinus arrhythmia. The presence of coupling 
between the cardiac and respiratory oscillatory processes 
has been revealed also using bispectral analysis in [5l|, [EH 
under both spontaneous and paced respiration. Our re- 
sults are in agreement with those when synchronization 
between the oscillating processes occurs as the result of 
their interaction. 



VI. CONCLUSION 

We have proposed the method for detecting synchro- 
nization from univariate data. The method allows one 



to detect the presence of synchronization of the self- 
sustained oscillator by external force with varying fre- 
quency. To implement the method one needs to analyze 
the difference between the oscillator instantaneous phases 
calculated at time moments shifted by a certain con- 
stant value with respect to each other. The instantaneous 
phases are defined at the oscillator basic frequency using 
continuous wavelet transform with the Morlet wavelet as 
the mother wavelet function. The necessary condition for 
the method application is the variation of the frequency 
of the driving signal. The method efficiency is illustrated 
using both numerical and experimental univariate data 
under sufficiently high levels of noise and inaccuracy of 
the basic time scale definition. 

We applied the proposed method to studying synchro- 
nization between the respiration and slow oscillations in 
blood pressure from univariate data in the form of R-R 
intervals. The presence of synchronization between these 
rhythmic processes is demonstrated within the wide time 
interval. The knowledge about synchronization between 
the rhythms of the cardiovascular system under paced 
respiration is useful for the diagnostics of its state (53|. 
The method allows one to detect the presence of synchro- 
nization from the analysis of the data of Holtcr monitor 
widely used in cardiology. 

The proposed method can be used for the analysis of 
synchronization even in the case when the law of the driv- 
ing frequency variation is unknown. If the frequency of 
the external driving varies in the wide range, the anal- 
ysis of the oscillator response to the unknown driving 
force allows one to make a conclusion about the pres- 
ence or absence of synchronization in the system under 
investigation. 



Acknowledgments 

We thank Dr. Svetlana Eremina for English lan- 
guage support. This work is supported by the Rus- 
sian Foundation for Basic Research, Grants 05-02-16273, 
07-02-00044, 07-02-00747 and 07-02-00589, and the 
President Program for support of the leading scien- 
tific schools in the Russian Federation, Grant No. SCH- 
4167.2006.2. A.E.H. acknowledges support from CRDF, 
Grant No. Y2-P-06-06. A.E.H. and A.A.K. thank the 
"Dynasty" Foundation for the financial support. 



[1] Blekhman I.I., Synchronization of dynamical systems 
(Moscow, Nauka, 1971). 

[2] Blekhman I.I., Synchronization in Science and Technol- 
ogy (American Society of Mechanical Engineers, 1988). 

[3] Pikovsky A., Rosenblum M., Kurths J., Synhronization: 
a universal concept in nonlinear sciences (Cambridge 
University Press, Cambridge, 2001). 

[4] Boccaletti S., Kurths J., Osipov C, Valladares D.L., 



Zhou C, Physics Reports 366, 1 (2002). 
[5] Rosenblum M., Pikovsky A., Kurths J., Schafer C, Tass 

P., in Handbook of Biological Physics (Elsiver Science, 

2001), pp. 279-321. 
[6] Meinecke F.C., Ziehe A., Kurths J., Miiller K.-R., Phys. 

Rev. Lett. 94, 084102 (2005). 
[7] Hramov A.E., Koronovskii A. A., Chaos 14, 603 (2004). 
[8] Hramov A.E., Koronovskii A. A., Kurovskaya M.K., 



10 



Moskalenko O.I., Phys. Rev. E 71, 056204 (2005). 
[9] Pecora L.M., Carroll T.L., Phys. Rev. Lett. 64, 821 
(1990). 

[10] Pecora L.M., Carroll T.L., Jonson G.A., Mar D.J., Chaos 
7, 520 (1997). 

[11] Pikovsky A., Rosenblum M., Kurths J., Int. J. Bifurca- 
tion and Chaos 10, 2291 (2000). 

[12] Boccaletti S., Pecora L.M., Pelaez A., Phys. Rev. E 63, 
066219 (2001). 

[13] Rulkov N.F., Sushchik M.M., Tsimring L.S., Abarbanel 

H. D.I., Phys. Rev. E 51, 980 (1995). 

[14] Pyragas K., Phys. Rev. E 54, R4508 (1996). 

[15] P. Tass, M. G. Rosenblum, J. Weule, J. Kurths, 

A. Pikovsky, J. Volkmann, A. Schnitzler, and H.-J. Fre- 

und, Phys. Rev. Lett. 81, 3291 (1998). 
[16] P. A. Tass, T. Fieseler, J. Dammers, K. Dolan, P. Mo- 

rosan, M. Majtanik, F. Boers, A. Muren, K. Zilles, and 

G. R. Fink, Phys. Rev. Lett. 90, 088101 (2003). 
[17] Chavez M., Adam C, Navarro, Boccaletti S., Martinerie 

J., Chaos 15 (2005). 
[18] Schafer C, Rosenblum M.G., Abel H.-H., Kurths J., 

Phys. Rev. E 60, 857 (1999). 
[19] Bracic-Lotric M., Stefanovska A., Physica A 283, 451 

(2000). 

[20] Rzeczinski S., Janson N.B., Balanov A.G., McClintock 
P.V.E., Phys. Rev. E 66, 051909 (2002). 

[21] Prokhorov M.D., Ponomarenko V.I., Gridnev V.L, Bo- 
drov M.B., Bespyatov A.B., Phys. Rev. E 68, 041913 

(2003) . 

[22] Hramov A.E., Koronovskii A. A., Ponomarenko V.L, 
Prokhorov M.D., Phys. Rev. E 73, 026208 (2006). 

[23] Hramov A.E., Koronovskii A.A., Levin Yu.I, JETP 127, 
886 (2005). 

[24] Hramov A.E., Koronovskii A.A., Physica D 206, 252 
(2005). 

[25] Hramov A.E., Koronovskii A. A., Popov P.V., Rempen 

I. S., Chaos 15, 013705 (2005). 

[26] Malpas S., Am. J. Physiol. Heart Circ. Physiol. 282, H6 
(2002). 

[27] Stefanovska A., Hozic M., Prog. Theor. Phys. Suppl. 139, 
270 (2000). 

[28] Adler R., Proc. IRE 35, 1415-1423 (1947). 

[29] Koronovskii A.A., Hramov A.E., JETP Lett. 79, 316 

(2004) . 

[30] Wavelets in Physics (Cambridge University Press, Cam- 
bridge, 2004), J.C. Van den Berg ed. 

[31] Koronovskii A. A., Hramov A.E., Continuous wavelet 
analysis and its applications (Moscow, Fizmatlit, 2003). 

[32] Lachaux J. P., Rodriguez E., Martinerie J., Varela F.J., 



Human Brain Mapping 8, 194, (1999). 
[33] Lachaux J. P., Rodriguez E., Van Quyen M.L., Lutz A., 

Martinerie J., Varela F.J., Int. J. Bifurcation and Chaos. 

10, 2429 (2000). 
[34] Le Van Quyen M., Foucher J., Lachaux J. P., Rodriguez 

E., Lutz A., Martinerie J., Varela F.J. J. Neuroscience 

Methods 111, 83 (2001). 
[35] Lachaux J. P., Lutz A., Rudrauf D., Cosmelli D., Le Van 

Quyen M., Martinerie J., Varela F. Neurophysiol. Clin., 

32, 157 (2002). 

[36] Quyen M. L., Foucher J., Lachaux J. -P., Rodriguez E., 
Lutz A., Martinerie J., Varela F. J. J. Neuroscience Meth- 
ods 111, 83 (2001). 

[37] Quiroga R.Q., Kraskov A., Kreuz T., Grassberger P. 
Phys. Rev. E 65, 041903 (2002). 

[38] Turalska M., Kolodziej W., Latka M., Latka D., West 
B.J. Shock 23, Suppl. 3, 90 (2005). 

[39] DeShazer D. J., Breban R., Ott E., and Roy R. Phys. 
Rev. Lett. 87, 044101 (2001). 

[40] Sosnovtseva O. V., Pavlov A. N., Mosekilde E., and 
Holstein-Rathlou N.-H. Phys. Rev. E 66, 061909 (2002). 

[41] Bandrivskyy A., Bernjak A., McClintock P. V. E., and 
Stefanovska A. Cardiovasc. Eng. Int. J. 4, 89 (2004). 

[42] Grossman A. and Morlet J., SIAM J. Math. Anal. 15, 
273 (1984). 

[43] Press W.H., Teukolsky S.A., Vetterling W.T., Flannery 
B.T., Numerical Recipes (Cambridge University Press, 
Cambridge, 1997). 

[44] Janson N. B., Balanov A. G., Anishchenko V. S., Mc- 
Clintock P. V. E. Phys. Rev. E 86, 1749 (2001). 

[45] Janson N. B., Balanov A. G., Anishchenko V. S., Mc- 
Clintock P. V. E. Phys. Rev. E 65, 036212 (2002). 

[46] Task Force of the ESC and NASPE, Circulation 93, 1043, 
(1996). 

[47] Shafcr C, Rosenblum M.G., Abel H.H. Nature 392, 239 
(1998). 

[48] Suder K., Drepper F.R., Schiek M., Abel H.H. Am. J. 
Physiol. 275, H33 (1998). 

[49] Hirsch J.A., Bishop B. Am. J. Physiol. 241, H620 (1981). 

[50] Kotani K., Hidaka I., Yamamoto Y., Ozono S. Method 
Inform Med. 39, 153 (2000). 

[51] Jamsek J., Stefanovska A., McClintock P.V.E., Khovanov 
LA. Phys. Rev. E 68, 016201 (2003). 

[52] Jamsek J., Stefanovska A., McClintock P.V.E. Phys. 
Med. Biol. 49, 4407 (2004). 

[53] N. Ancona, R. Maestri, D. Marinazzo, L. Nitti, M. Pel- 
licoro, G.D. Pinna, S. Stramaglia, Physiological Measure- 
ment 26, 363-372 (2005). 



